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ABSTRACT 

One of the more significant error sources in using the 
for attitude determination of aircraft and earth orbiting satellites is multipath interference. 
The GPS signal is reflected from the spacecraft structures resulting in a composite signal 
that differs from the direct line-of-sight signal in terms of phase and amplitude. The 
resulting phase errors can have a significant impact on the attitude determination error. 
This thesis will describe the GPS signal generation and provide tools to simulate the 
signal. A Kalman filter was developed to track the phase changes in the simulated GPS 
signals and its performance described in the absence and presence of multipath. 
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I. 



INTRODUCTION 



A. BACKGROUND 

In order for a spacecraft to perform its mission, it must be able to orient itself in 
relation to either the sun, the earth or both within an inertial reference frame. The primary 
system used to accomplish this is the spacecraft attitude determination and control 
subsystem (ADCS). The ADCS uses external references to determine the vehicle’s 
absolute attitude with respect to some fixed inertial reference frame. Traditionally, these 
external references have included the Sun, the Earth’s infrared (ER) horizon, the local 
magnetic field direction, and the stars. Recent technological developments have added a 
new tool to the external reference arsenal - the Global Positioning System (GPS). 

Typical traditional ADCS sensor characteristics are shown below in Table 1. 
(Larson and Wertz) 



SENSOR 


TYPICAL PERFORMANCE 
RANGE 


WEIGHT (KG) 


POWER (W) 


Inertial measurement unit 
(gyros & accelerometers) 


Gyro drift rate = 0.003°/hr to 
l°/hr 


3 to 25 


10 to 200 


Sun Sensors 


Accuracy=0.005° to 3° 


0.5 to 2 


0 to 3 


Star Sensors (scanners & 
mappers) 


0.0003° to 0.01° 


3 to 7 


5 to 20 


Horizon Sensors 


O 

r— l 

2 

O 

O 


2 to 5 


5 to 10 


Magnetometer 


0.5° to 3° 


0.6 to 1.2 


< 1 



Table 1 Typical ADCS Sensors 



In comparison, the Trimble TANS Vector attitude determination system using 
GPS is capable of an accuracy of 0. 1° with a mass of 2.2 kg (GPS receiver and antennas) 
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and a power requirement of 7.5 W. Additionally, the Trimble system also provides three- 
dimensional orbital position information; which other sensors are incapable of providing. 

As described by Cohen (1996), attitude determination using GPS is achieved by 
determining the relative positions of multiple antennas, mounted to the vehicle, based on 
differences in the carrier phase measurements at the different antennas. The basic premise 
is that the GPS satellite is so distant, relative to the receiving antenna separations, that 
arriving wavefronts can be considered as effectively planar as shown below in Figure 1 . 



The phase of a signal, <f > , is a function of the wavelength, X, and distance traveled, 



The wavelength is a function of the speed of light, c, and the frequency, f. The 
distance traveled is a function of the speed of light and the time difference of arrival 
between antennas. At: 




Figure 1 Attitude Geometry 



d: 
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d = c At 



Therefore a relationship between, time and phase can be shown to be 
C At 




The incoming signal will arrive at the antenna closest to the transmitter before 
reaching the most distant antenna. By measuring the difference in carrier phase between 
the antennas, a receiver can determine the relative range between any two antennas. The 
measured differential phase, A <p , for baseline /' and satellite /, is proportional to the dot 
product of the baseline vector, x , and the line of sight unit vector to the transmitter, s . 
As shown in Figure 2, the GPS receiver measures only the fractional part of the 
differential phase. The integer component, k, must be resolved by independent means. 

The differential phase can therefore be expressed as: 



measurements between the different antennas, the spacecraft attitude can then be resolved 
as shown by Cohen (1995). Filtering this measurement will be the focus of this thesis. 

The dominant noise source in differential phase measurements is multipath, which 
occurs when the signal arriving at the antenna consists of the reflected signals in addition 



transmission and appears as additive noise at the antenna. Because the antenna locations 
are different, the multipath signature at each antenna is unique and will not be common to 
every antenna. Multipath accounts for more than 90% of the total error budget in carrier 
phase measurements. (Lightsey) 




where Vy is the noise in the measurement. Given the differential phase 



to the line-of-sight source. The reflected signal is phase shifted with respect to the original 
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B. PURPOSE 



This thesis will discuss the generation of the GPS carrier signal which is used in 
making differential phase measurements. Models will be developed for simulating the 
multipath environment and the GPS signal. A Kalman filter which tracks the phase in the 
signal for use in the differential phase measurements will be described. Finally, results of 
the Kalman filter’s use in a multipath environment will be presented. 

Chapter II will discuss the GPS Signal Generation and present a model of the 
signal. The multipath interference will be examined and modeled in Chapter III. Chapter 
IV provides the derivation of the Kalman filter for tracking the phase of a GPS signal. 
Chapter V will show the results when the filter is used in tracking the GPS phase. The 
final chapter will summarize the findings and present other ideas for minimizing multipath 
and suggestions for future investigation. 
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H. GPS SIGNAL GENERATION 



The GPS signal is generated as shown in Figure 3. The GPS satellite transmits 
two carrier frequencies, LI (primary frequency) and L2 (secondary frequency). These 
carrier frequencies are modulated by spread spectrum codes with a pseudorandom noise 
(PRN) sequence unique to each satellite. The LI frequency is modulated by the 
navigation message along with two PRN codes: the coarse/acquisition (C/A) code and the 
precision (P) code. The GPS L2 signal is modulated by one of the codes using Binary 
Phase Shift Keying (BPSK). The GPS LI signal is modulated using both codes by an 
unbalanced quadrature phase shift key (QPSK) modulator which is actually two BPSK 
modulators with different amplitudes. This development closely follows that given by 
Kaplan (1996). 

A. PN CODE GENERATION 

The PRN codes used in the generation of the GPS signal are non-maximum length 
codes also known as Gold codes. These codes are generated from an n-bit shift register 
generator. This register has an initial state or fill. After each clock cycle the output of a 
specified cell is used as the input to a modulo-2 adder (exclusive-or). The output of the 
modulo-2 adder is then fed back into a specified stage of the register. As an example, 
consider the shift register shown in Figure 4. This register uses the 3 rd and 5 th stages as 
the input to the modulo-2 adder and feeds the result back into the first stage. The 
polynomial used to describe the design of linear code generators is of the form 

i+2>' 

where x* means the output of the i* cell of the shift register is used as the input of 
the modulo-2 adder and the 1 means the output of the adder is feed in to stage 1 . The 
polynomial for the register in Figure 4 is then 1 + X 3 + X 5 . If this register was initially 
filled with all 1 ’s at time 0, the resulting register for the first ten clock cycles would be as 
shown in Table 2. The binary signal output is taken from the last register and in this case 
would be 1 1 1 1 1000110. 
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Information 



Figure 3 GPS LI and L2 Signal Generation 



Clock 




Output 



Shift Register 

Figure 4 PRN Code Shift Register for Example Polynomial 1+xr+x 5 
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Register 

time 1 2 3 4 5 

0 11111 

1 0 1111 

2 0 0 1 1 1 

3 0 0 0 1 1 

4 1 0 0 0 1 

5 110 0 0 

6 0 110 0 

7 10 110 

8 110 11 

9 1110 1 

10 0 1110 

Table 2 PRN Code Generation Example Results for Polynomial l+x 3 +x J 



The Matlab m-function gen poly.m can be used to generate the PRN sequence for 
any user defined polynomial and register. 

B. C/A CODE GENERATION 

The GPS coarse/acquisition (C/A)-code is a Gold code generated by two 10-bit 
shift registers, G1 and G2. The polynomial for each is given below (Kaplan, 1996) 

Gl = 1 + X 3 + X 10 

G2 = 1 + X 2 + X 3 + X 6 + X 8 + X 9 + X 10 

Both C/A-code registers have an initial fill of all 1 ’s. Each code is derived and the 
result of G2 is then delayed by the satellite PRN code number. This delay effect is 
obtained by the exclusive-or of the selected positions of the G2 register. The C/A-code is 
then obtained by the exclusive-or of the Gl and delayed G2 signals. An example for 
satellite vehicle 1 is shown in Figure 5. 

The code tap positions and initial code sequences are given in Table 3. The m-files 
gl. m and g2.m generate the codes for a specified satellite vehicle number. 
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C/A Code 
Gi(t) 



Figure 5 C/A-Code Generator for Satellite Vehicle 1 with C/A Code Tap Selection 2 @6 

(from Kaplan) 
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1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 



C/A-Code 

Tap 

Selection 


C/A-Code 

Delay 

(Chips) 


P-Code 

Delay 

(Chips) 


First 1 0 C/A 

Chips 

(Octal) 


First 12 P- 

Chips 

(Octal) 


2©6 


5 


1 


1440 


4444 


307 


6 


2 


1620 


4000 


408 


7 


3 


1710 


4222 


509 


8 


4 


1744 


4333 


109 


17 


5 


1133 


4377 


2010 


18 


6 


1455 


4355 


108 


139 


7 


1131 


4344 


209 


140 


8 


1454 


4340 


3010 


141 


9 


1626 


4342 


203 


251 


10 


1504 


4343 


304 


252 


11 


1642 


4343 


506 


254 


12 


1750 


4343 


607 


255 


13 


1764 


4343 


708 


256 


14 


1772 


4343 


809 


257 


15 


1775 


4343 


9010 


258 


16 


1776 


4343 


104 


469 


17 


1156 


4343 


205 


470 


18 


1467 


4343 


306 


471 


19 


1633 


4343 


407 


472 


20 


1715 


4343 


508 


473 


21 


1746 


4343 


609 


474 


22 


1764 


4343 


103 


509 


23 


1063 


4343 


406 


512 


24 


1706 


4343 


507 


513 


25 


1743 


4343 


608 


514 


26 


1761 


4343 


709 


515 


27 


1770 


4343 


8010 


516 


28 


1774 


4343 


106 


859 


29 


1127 


4343 


307 


860 


30 


1453 


4343 


308 


861 


31 


1625 


4343 


409 


862 


32 


1712 


4343 



Table 3 GPS Code Phase Assignments (from GPS-ICD-2000) 
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test_CA.m is a script file to test C/A code generation. This program calculates 
the first 10 C/A-chips and converts the result to octal numbers that match the First 10 
C/A-Chips (Octal) of the GPS-ICD-2000 table shown in Table 3 

For example the first 10 chips of the PRN 5 C/A-code is 100101 1011 which is 
1133 in octal. 

C. P-CODE GENERATION 



The GPS Precision (P) Code is generated by PRN sequences using four 12 bit shift 
registers whose polynomials and initial states are shown below: 



Register 


Polynomial 


Initial State 


X1A 


1 +X 1 +X 2 +X 5 +X 8 + X 9 +X 10 +X 11 + X 12 


001001001000 


X1B 


1 +x 1 +x 2 +x 3 +x 4 +x 5 + x 7 +x 8 +x 9 +x 10 +x 11 +x 12 


010101010100 


X2A 


1 + X* + X 3 + x 4 -!- X 5 + X 7 + X 8 + X 9 + X 10 + x 11 + X 12 


100100100101 


X2B 


1 + X 2 + X 3 + X 4 + X 8 + X 9 + X 12 


010101010100 



Table 4 P-Code Generator Polynomials and Initial States 



The P-code is 7 days in length at a chipping rate of 10.23 MBps. The 7 day 
sequence is the Modulo-2 sum of two subsequences XI and X2i; their lengths are 
15,345,000 chips and 15,345,037 chips respectively. 

XI itself is generated by the Modulo-2 sum of the output of two 12-stage registers 
(XI A and X1B) short cycled to 4092 and 4093 chips respectively. When the XI A short 
cycles are counted to 3750 the XI epoch is generated. This occurs every 1.5 seconds, 
after 15,345,000 chips of the XI pattern have been generated. The XI period is defined 
as 3750 XI A cycles (15,345,000 chips) which is not an integer number of X1B cycles. To 
accommodate this the X1B shift register is held in the final state (chip 4093) of its 3749 th 
cycle for 343 chips. The X2 sequence is generated in a similar manner, however the X2 
period is defined to be 37 chips longer than the XI period in order to force the X2A and 
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X2B epochs to precess with respect to the XI a and X1B epochs. This precession 
continues until the end of the 7 day period. 

The Y-code 1 is used when the anti-spoofing mode of operation is activated. It is 
generated in the same fashion as the P-code, however the polynomial generating equations 
are classified. (ICD-GPS-2000) 

xla.m, xlb.m, x2a.m and x2b.m are used to calculate the PN codes. The m-file 
test_P.m is used to test P-code generation. 

D. DATA GENERATOR 

The output of the data generator is the navigation message, D(t ) . The message 
includes satellite vehicle ephemerides, system time, clock bias information, status 
messages, and C/A to P-code handover information. The navigation message is generated 
at 50 bps and is added to the P and C/A-codes to modulate both the LI and L2 signals. 
(ICD-GPS-200) 

E. BPSK MODULATOR 

After the various codes are generated they are then used to modulate the carrier 
signal using binary phase shift keying, BPSK, resulting in the LI or L2 signal. The 
Simulink program used to simulate the BPSK modulator is shown below in Figure 6. 



1 The Y-code and P-code are sometimes discussed together and are referred to as P(Y)-code meaning P (or 
Y) code. 
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Inport 



■►JV — HU- 



Zero Order Hold 



Relay 



B 



Sine Wave 



Product 



Mux 



BPSK 

To Workspace 



Mux 



Outport 



Figure 6 BPSK Simulink Diagram 

The PRN code is input via the inport. It is sent through a zero-order-hold and 
relay function to generate a non-return to zero (NRZ) signal, which is multiplied by the 
sine wave, resulting in a binary phase shifted signal. 



14 



BPSK Simulation 




time 



Figure 7 BPSK Simulation Results 



Figure 7 shows an example of the BPSK simulation results. This figure shows the 
results: NRZ signal in the top graph, the sine wave in the second graph, the product of the 
two, or the BPSK signal in the third graph, and a composite of the data and resulting 
signal in the bottom graph. 



F. LI SIGNAL 

The P-code at, 10.23 x 10 6 chips per second, and C/A-code, at 1.023 x 10 6 chips 
per second, are combined with the 50 bits per second data and are then used to modulate 
the LI frequency ( 1 54 • f 0 = 1 575.42 MHz ) as shown in Figure 3 . This results in an 
Unbalanced Quadri-Phase Shift Keying (UQPSK) modulation with the data bits added to 
the ranging codes, the C/A-code signal lagging the P-code signal by 90°; and the C/A- 
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code signal power nominally exceeding the P-code signal power by 3 dB (Struza) 
Mathematically, this can be represented by: 

£,(<y,0 = 4[-P(0® £>(0] c os(fi) 1 0 + V2i4,[G(/)®i)(/)]sin(<o 1 /) 
where P{t) is the P-code, D(t) is the data, and G(t ) is the C/A-code. 

G. L2 SIGNAL 

The L2 frequency (120- f a = \221.6MHz) is modulated by either the exclusive-or 
of the P-code and data, the P-code, or the exclusive-or of the C/A code and data, as 
selected by the control segment. Mathematically, this can be represented by: 

L 2 (a> 2 t) = ^ [/*(*) ® D(0]cos(<y 2 0 
L 2 (co 2 t) = A 2 [P{t)]cos(co 2 t) 

L 2 {co 2 t) = ^ 2 [G(0]cos(tu 2 /) 



H. SUMMARY 

The GPS signal consists of two separate signals, LI and L2. These signals are 
generated using BPSK and unbalanced QPSK modulation. The modulation signal is 
generated using the P and C/A codes to which the data signal is added. The Matlab files 
init.m, sig gen.m and my glb.m are used to initialize and generate a sample of the codes 
used to modulate the signal. These codes are then used by the Simulink programs 
LI sim. m and L2 sim.m to generate a sample of the signal. In the next chapter the 
multipath environment will be described. 
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m. MULTIPATH 



A. INTRODUCTION 

This chapter will discuss the multipath environment that impacts the use of GPS 
for attitude determination and will present a method to model multipath interference using 
Matlab. 

Multipath occurs when the signal arrives at the antenna from reflected surfaces and 
from the line of sight sources. A simple case with a single reflected path signal is shown in 
Figure 8. 




The reflected signal is phase shifted with respect to the original transmission and 
appears as additive noise at the antenna. When antenna locations are different, the 
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multipath signature at each antenna is unique and the error in the phase measurement is 
not common. (Lightsey) Multipath interference, due to environmental factors, will impact 
the signal, however this effect will be common to each received signal. 

For the developed model it will be assumed that the multipath environment will be 
static. In reality the multipath environment will be varying due to the satellites attitude 
and position relative to the GPS satellites changing over time. 

B. MATHEMATICAL MODEL 

In the absence of the multipath, the input signal to a GPS receiver for the LI 
QPSK signal is given by: 

s(t)= A c cos(cot +a(t)f) 

where A* is the received signal amplitude and a(t) = 0,1,2 or 3 is the pseudo- 
random code waveform that phase modulates the carrier. (Kumar) In the presence of N 
multipaths, in addition to the direct line of sight path, the input signal is given by 

s m (t) = a 0 s(t ) + a } s(t - r, )+. . .+a w s(? - r N ) 

where a i and r i denote the amplitude and delay of the i* multipath for 
i= 1 ,2, . . ,,N. Graphically, this can be represented as shown in Figure 9. 
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cc N 



s m (t) 



Figure 9 Multipath Block Diagram 
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Typically for BPSK and QPSK signals the carrier frequency operates at a much 
higher frequency than the modulation rate. As a result, a small value of delay, t , will 
result in large phase differences. When the multipath signal with a large phase difference 
is added to the direct path there may be destructive adding which can significantly impact 
the received signal. 

C. SIMULATION MODEL 

In the model mpath.m, three vectors are used to develop the signal s m 

delay -1 by z N vector of delayed samples. The first entry is the undelayed 
signal, the next entry is the signal delayed by one sample period. The length of this vector 
is the maximum value in the tau vector and is initialized as all zeros. 

tau - vector of delay times. Each entry is 1 minus the delay time in sample 
periods, dt, and is used to index the delay vector. The first entry is 1 and represents a 
delay of zero. 

alpha - vector of amplitudes. Each entry is the amplitude of the 
corresponding delayed signal. The first entry is a 1 in order to return the undelayed fine 
of sight signal. 

As an example, the BPSK signal shown in Figure 10, with N=3 multipath signals, 
will be used to demonstrate the use of the vectors and the algorithm used to develop the 
signal 

S m ( 0 = «o s (t - T 0 ) + cc x s{t -r 1 ) + a 2 s(t-r 2 ) + a 3 s(t - r 3 ) 

where 

r = [l 10 12 20] 

a = [ 1 0.5 0.3 0.1] 

thus 

s m (t) = s(/) + 0.5*s(f-9**ft) + 0.3*s(f-ll‘<#) + 0.1*s(/- \9-dt) 
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If the sample period, dt, equals 0.01 the following signal is produced using the m- 
file mpath.m 



Multipath simulation 
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Figure 10 Multipath Example Signal 

The top graph in Figure 10 shows the undelayed signal, the middle graph shows 
the direct path and multipath signals plotted individually and the bottom graph shows the 
sum of the individual and multipath signals. 

If the multipath delay time and amplitude are increased with the following 
parameters 
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combined composite undelayed 



r = [l 20 22 30] 

a = [l 0.5 0.3 05] 

the resulting multipath signal will be as shown here: 
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Figure 11 Multipath Example with Large Delay 

Notice that in both examples the multipath signal has a different amplitude from 
the original signal and that the larger the delay the larger the change in the phase of the 
signal Below is an example of destructive interference obtained with the values 
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combined composite undelayed 



r = [l 50 52] 
a = [l 0.5 0.4] 
Multipath simulation 




Figure 12 Multipath Example with Destructive Interference 

This chapter has described the multipath environment and presented a simulation 
model to generate a signal based on user defined parameters. The next chapter will 
describe a Kalman filter which is used for tracking the GPS signal. 
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IV. KALMAN FILTER FOR TRACKING PHASE 



A. DISCRETE KALMAN ESTIMATOR 

The development of a filter to track the GPS signal in a multipath environment 
began with the implementation of a Discrete Kalman Estimator for tracking the phase 
changes of BPSK signal. The derivation below closely follows that given by Professor 
Kirk (1975). 



1. Plant Equations 



The plant equations were derived as follows: 
Xj = cos(tu t) 



x 2 = x, = -co sin(ty t ) 



x(t) = 






\x 2 J 



r cos {cot) ^ 
V-<y sin(<y t)J 



i(0 = 






\x 2 J 



r -<ysin(dj/) ^ 




f o r 


r cos {cot) ^ 


V-6J 2 cos {cot)j 




K-CO 1 Oo 


V-w sin(o; t)J 



= Ax(0 



2. Plant Characterization 



The discrete model of the plant is characterized by the linear discrete difference 
equations: 



x(k + 1) = tf> x(k) + w(k) 



z(k) = Cx(k) + v(k) 

where: 

x(k) is the n-dimensional state vector at time k, 

z(k) is a q-dimensional output measurement vector at time k 

w(k) is a p-dimensional vector of random forcing inputs at time k 
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<}> is the state transition matrix 
C is the observation matrix 
and 

v(k) is the q-dimensional vector of random measurement noise at time k 

3. Gain, and Covariance Equations 

The gain, G, and covariance equations derived by Kirk (1975) are: 

G(k) = P(k\k - 1) • C T • [C • P(k\ k-\)C T + /?]"' 

P{k\k) = [I - G(k) ■ C] ■ P(k\k - 1) 

P(k + \\k) = </>-P(k\k)-0 T +O 

where 

C is the observation matrix 
I is the identity matrix, 

G is the Kalman gain 
P is the variance of estimation error 
Q is random process forcing input covariance matrix 
R is the variance of measurement noise. 

The estimator equations are: x(k\k) = x(k\k - \) + G(k)-[z(k)~ C • x(k\k - \)\ 2 
x(k\k - 1 ) = <f>- x(k - 1 | k - 1) + A -u(k - 1) 

with the initial condition jc(0| — 1) = x o 



2 * s rea d as x hat of j given j and is defined as the estimate of x at time j given measurements up to and 

including time j. 
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4. Assumptions 



a) The measurement noise has zero mean. 

£[v(i)]=0 k = 0,1,... 3 

is uncorrelated and has covariance R(k): 

£[v(*)v r O')] = £[vO> r (*)] = R(k)S, 
where S jk is the Kronecker delta function defined by 

o _Jl J = * 

Jk [0 j*k 

b) The initial state is a random variable with known mean 

4<o)]=r 

and covariance: 

i{[*(0) - * 0 ] [x(0) - x a ] T j = M 

The measurement noise and initial state are uncorrelated: 

£[x(0)v r (£)] = £[v(*)x r (0)] =0 k = 0,1, . . . 

c) by the linear relationship 

x(k\k) = x(k\k -\)+G(k)- [z(k) -C-x(k\k- 1)] 

where G(k) is the Kalman gain and z(k) - C ■ x(k\k - 1) is the correction term, 

the difference between the measured and predicted value. 

d) The random process forcing input has zero mean 

£[w(£)] = 0 k = 0 , 1 ,... 

is uncorrelated and has covariance Q(k) 

E[»(k)w T to] = E[m,(J> T W] = QWS, j,k = 0,1,2,... 

e) The forcing random process and initial state value are uncorrelated 

3 £[] is the expected value or mean of the random variable. 
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£|w(£);c r (0)j = £|x(0)w r (£)j = 0 k = 0,1,2,... 
f) The forcing random process and measurement noise process are uncorrelated 
£[w(i)v r 0)] = £[v(y)w r (t)] = 0 k = 0,1,2,... 

COMPUTATIONAL PROCEDURE 

The computational procedure is shown in the following algorithm. 

Initialize with 
k=0 

x(o|-i) = 

P(k\k-l) = M 

The procedure iterates on k in the loop: 

G(k) = P(k\k - 1) • C T • [C ■ P(k\k - 1) • C T + tf]"' 

P(k\k) = [I - G(k ) • C]P(k\k - 1) 

P(k + 1|£) = ip- P(k\k) -<p T +Q {this will be P(k\k - 1) next time through) 
take measurement z(t) 

x{k\k) = x(k\k - 1) + G(k) ■ [z(k) - C ■ x(A:| k - 1)] 

£(£|A:-1) = <f>-x{k -\\k -\) + k-u{k -X) 
k=k+l 

End of the Loop 



PHASE TERMS 



The QPSK signal can be represented by 



s(t ) = A cos (cot + <j> (()) 

where <f>(t) represents the phase change of the signal 



t(0 = 



to 

tc+- 
y ° 2 

to+X 

to + 

2 



Substituting in these phase values and simplifying using trigonometric identities 
results in 




Acos(a>t + <f> 0 ) 
-As\n{(ot + <(> 0 ) 
-Acos(cot +<j> 0 ) 
Asm(cot + <f> c ) 



Using these relationships and the sine and cosine information in the state estimates, 
the Kalman filter can estimate which phase change has occurred. The error term is 
calculated by finding the estimate of the signal that minimizes the difference between the 
measured value. Assuming 

r cos (cot + <f> 0 ) ^ 



x(k\k - 1 ) = 



-cosin(cot + 0 o )J 



C = [1,0] 

C p , = [0,l] 



the four error terms can be calculated as shown below, and the minimum value is 
used to update the prediction 

e 0 = z(k ) - Cx(£|/r - 1) = z(k) - cos(*y k + <j> 0 ) 
e } = z{k ) + C ps \{k\k — 1) / co — z(k ) - sin(<y k + <j> o ) 
e 2 = z(k) + Ci(k\k - 1) = z(k) + cos(a> k + <f> 0 ) 
e % = z(k)~ C ps \(k\k - 1) / at = z(k) + sin(<u k + <j > 0 ) 
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These error terms are then used to update the estimate of the signal without any 
phase shifting which will be called y, and an estimate of the phase shifted signal called 
y_ps. 

Next, the results of using the developed Kalman filter to track the phase of a 
BPSK, QPSK, and unbalanced QPSK signal will be shown. 
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V. 



INITIAL RESULTS 



This chapter will show the results obtained when the Kalman filter, discussed in the 
previous chapter, is used the track the phase of a BPSK, QPSK and unbalanced QPSK 4 
signal. The performance of the Kalman filter, will be shown with and without noise, and 
in the presence of multipath interference. It will be shown that the Kalman filter will track 
the combined line of sight and direct path signals phase and amplitude. It will then be 
shown that when given a priori information the filter can distinguish between the line of 
sight and multipath signals and correctly track them, given a static multipath environment. 

A, SIGNAL GENERATION 

The Matlab m-file uivars.m sets up the user interface shown in Figure 13 that can 
be used to control the simulation. 

Each entry can be edited by the user to change each of the parameters discussed 

below: 

• Signal freq - Sets the frequency, of the oscillator, in Hertz, used by the 
carriers in the BPSK and QPSK signal generation. 

• Phase Offset - Determines the initial phase offset, of the line of sight signal in 
degrees. 

• Amplitude - Sets the amplitude of the generated signal. 

• Ts - This sets the sampling period in seconds. In the example below Ts =0.01, 
resulting in 100 samples per second. 

• Stop time - Sets the Simulink parameter stop time to control length of 
simulation. 

• Tb - Sets the bit period which determines the length of time for zero-order hold 
per data bit. 



4 An unbalanced QPSK signal is a QPSK signal with differing energy in the in-phase and quadrature 

phase components. For GPS the LI signal the amplitude of the quadrature term is equal to V2 times the 
in phase term. 
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SNR dB - Using this value and the entered amplitude value of the signal, the 
random number generator is multiplied by a scaling factor to simulate the 
desired SNR level. 

tau - Delay vector used to control the multipath delay times. See Chapter III 
Section C for more information. 

alpha - Amplitude vector - Used to control the multipath amplitudes. 




Figure 13 User Interface to Control BPSK and OPSK Simulation 




The pulldown menu Init initializes all the entered values. Plot gives the user a 
choice of plots to choose from to document the output of the simulation. Sim allows the 
user to choose which simulation will be run. 

B. BPSK 

Using Matlab and Simulink, a BPSK signal was generated and the algorithm 
discussed in the previous chapter was used to track the phase changes in the signal. 5 

1. BPSK Signal Without Noise 

Figure 14 shows the performance of the algorithm when there is no noise present. 
In this graph, x is the plant state, z is the measured signal (x +noise), y is the recovered 
carrier, and y_ps is the recovered phase shifted signal. The last graph shows the error in 
the recovered signal which is the difference in the generated signal and recovered estimate, 
x-y_ps. Notice there is an initial error during initialization, however the error quickly goes 
to zero as expected in the absence of noise. 

2. With Noise 

Keeping all variables the same as above and adding zero mean Gaussian white 
noise, by using the random function in Matlab, with an SNR value of 12 dB 6 , the carrier is 
still recovered, however small errors are now introduced in the estimate of the phase 
shifted signal. 



5 See Matlab m-files kaljnit.m, kalman2.m and the Simulink program kaljsim 

6 Spilker states that typical values for C/N 0 are 40.6 dB-Hz, which corresponds to an SNR of 25 dB. See 
Appendix C for a link budget estimate of the signal to noise ratio. 
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Figure 14 Kalman Filter as Phase Lock Loop for BPSK Without Noise 
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Kalman Filter 




time ~ sec 

Figure 15 Kalman Filter for BPSK with Noise 



C. QPSK 

Next the results for a QPSK simulation are shown. These plots were generated 
using kalqsim.m 

1. Without Noise 

Figure 16 shows the performance of the algorithm when there is no noise present. 
The true signal, x, and the estimate y_ps and the error in the estimate, x-y_ps are all 
plotted. 
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Kalman Filter Error 
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Figure 16 Kalman Filter for QPSK, Without Noise 



2. With Noise 

Figure 17 shows the estimation when white noise 7 is added with an SNR value of 
12 dB. Figure 18 shows the error in the estimation. The carrier is still recovered, however 
small errors are now introduced in the estimate of the phase shifted signal. The cause of 
the error will be explained in the next section. 



Throughout this study only gaussian white noise was considered. Colored noise from other sources, such 
as electromagnetic interference was not considered. It was assumed that any such noise would have the 
same affect at both receivers and would not change the difference in the measured phase values. 
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Kalman Filter 




Figure 1 7 Kalmcm Filter for OPSK with Noise 




Figure 18 Error in OPSK Estimator 
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D. SOURCES OF ERROR IN NOISY CASE 



The error in the estimation caused by the noise can best be explained by examining 
the BPSK simulation. In the algorithm, the Kalman filter tracks the carrier signal and 
estimates which generated phase shifted version has the minimum error when compared 
with the measured signal z. Figure 19 is a plot of the original BPSK signal, x, the signal 
with noise added, z; the estimated signal, y_ps; the two possible signals, xO and x2, that 
must be decided between; and the error, eO and e2, used as the basis for the decision. 

The error in the signal is best described as “clicks” which occur when the Kalman 
filter is unable to distinguish between the correct and incorrect phase because they both 
have approximately the same value. The reader is invited to examine this effect between 
times 1 and 2 seconds. The correct estimated signal is X2 which results in the minimum 
error over the bit period, however at approximately time 1.9 the out of phase signal 
estimate XO has the same value as X2 and the error eO is less than el due to the added 
noise. This produces an erroneous decision to switch phase. Typically the erroneous 
decisions all occur when both the correct and incorrect signal estimates are near a zero 
crossing. This is due to the nature of the cosine function which results in 
cos(90°)=cos(270°)=0. 
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Kalman Filter - Signal Options 
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Figure 19 Kalman Filter Error Terms and Decision Options 
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E. UNBALANCED QPSK 



As discussed in Chapter II, the GPS LI signal is actually an unbalanced QPSK 
signal with differing energy in the quadrature and in-phase components. As shown in 
Chapter II Section F on page 15 the LI signal is given by 

L x (<y, t) = Aj [P(t) ® D(/)] cos(<yjQ + *JlA l [G(Q ® £>(/)] sin(<y,Q 
An example of the unbalanced QPSK signal is shown in Figure 20. The top graph 
shows the composite of the binary data representing [G(Q ® D(t)] with the quadrature 

component, the middle plot shows the binary data representing [P^)® D(t)] with the in- 

phase component. In both cases the Ai has been normalized to 1 . The last plot shows the 
sum of the two signals which is the resulting unbalanced QPSK signal. 



QPSK Simulation 




Figure 20 Unbalanced QPSK Signal Generation 
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When the Kalman filter is used to track this signal as shown in Figure 21, errors 
occur when the in-phase and quadrature terms have differing data. The filter locks on to 
the initial amplitude and assumes the amplitude will not change. In this example the 
absolute value of the amplitude is the same when the data bits are the same (both +1 or 
both -1) and is different with the data bit have opposite values. The state estimation was 
developed under the assumption the signal would have a constant amplitude. In order to 
resolve this problem the logic for phase change detection would need to be modified. 
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Figure 21 UOPSK Without Noise 
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When noise is added to the unbalanced QPSK as shown in Figure 22 the 
performance is the same as in the QPSK example when the data bit for the in-phase and 
quadrature phase term is the same, however when they are different, the performance 
degrades. The error in the estimation is shown in Figure 23. 
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Figure 22 UQPSK with Noise 
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Kalman Filter Error 





Figure 23 Error in UQPSK Estimator 



F. MULTIPATH RESULTS 

1. Without A Priori Knowledge 

The multipath parameters used to generate Figure 10 are used to generate the 
measured signal z and the results obtained are shown below in Figure 24. The estimate, 
y_ps, and the direct path signal, x, are shown plotted together and the error, or the 
difference between x and y_ps is also shown. 
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Kalman Filter Error 




time ~ sec 



Figure 24 Results of Kalman Filter Estimation in Presence of Multipath 



The Kalman filter has correctly determined the phase changes, however due to the 
multipath, the amplitude of the estimation is now different from the original signal. Using 
the larger multipath spread parameters used to generate Figure 1 1 the resulting phase 
error is observed as shown Figure 25. 
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Kalman Filter Error 




Figure 25 Results of Kalman Filter Estimation in Presence of Multipath with Increased 

Delay 



The results of these last two plots indicate the Kalman filter performance in the 
presence of multipath is much worse than in the presence of white noise. If multipath is 
present it will introduce a phase error that will be unique to each antenna. When the 
differential phase is calculated it will now include the multipath phase error and will no 
longer be just a function of the geometry. 
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2 . 



With A Priori Knowledge 



If the multipath amplitude and time delay are known, the Kalman filter can be 
modified to track both the the direct path signal and delayed multipath signals. A simple 
case with a single multipath channel is shown below with the following assumptions: 

The combined signal is given by: 

s m( t ) = s(t) + as(t-r) 

The incoming signal is not changing phase and is given by 

s(t ) = A cos(<y t + 9(t)) 

where 0(t)= a(t) f and a(t) = 0,1, 2, or 3 representing the intelligence for the BPSK 
or QPSK modulation. 

The time delay r is converted into the phase delay <j> 

s(t-r) = A cos(<y(/ - r) + 9(t - r)) = A cos(<y t - co r + 9(t - r)) = A cos(<y t + <f>) 
and the composite multipath signal is now given by 

s m ( t ) = A cos(<u t + 9) + a A cos(<y t + ^ ) 

The plant equations are now 

Xj - A cos {co t + 9) 

x 2 = x, = -co A sin(«y t + 9) 

x 3 = a A cos (co t + <j>) 

x 4 = x 3 = -co a A sin(<y t + <t > ) 



x(/) = 







r A cos(<y t + 9) ^ 


x 2 




— co A sin(tu t + 9) 


*3 


— 


a A cos(<y t + <j>) 






co a A sin(<y t + cf>)) 



f*0 




f -co A sin(n) t + 6) ^ 




"0100^ 




f A cos(<y t + 9) ^ 


*2 




-co 2 Acos(cot9) 




-co 2 0 0 0 




- coA sin(<u t + 9) 


*3 




- co a A sin(fi> t + </> ) 




0 0 0 1 




a A cos (co t + <f>) 


U J 




k- co 2 a A cos(co t + <f>)j 




v 0 0 -co 2 Oy 




co a As\n(co t + <f>)) 



= Ax(/) 
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Since a priori knowledge is known the values for the estimate at time zero can 
now be initialized as follows: 



x(0|-l) = 



A 

0 

i a. Acos(<f>) 
-aa A sin(^) 



Figure 26 shows the line of sight signal, delayed signal, and the measured signal z 
which were generated using a=0.5 and t=10 with Ts=. 1. For this example the multipath 
environment is assumed to be static. In reality, the multipath environment would be much 
more dynamic due to the satellites orbit, the satellites changing attitude, and the changing 
position of the GPS satellites. 





N 




Figure 26 Example Signal for A Priori Case 
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los & xkkps(1 ,:) (estimate = 




Figure 27 A Priori Knowledge Results Showing Signals and Estimates 



Figure 27 shows the estimations and original signals plotted on top of each other. 
This demonstrates how, given a priori knowledge, the line of sight and delayed signal can 
be tracked. Now there is no phase error due to the multipath signal, unlike the examples 
shown without a priori knowledge when the filter tracked the phase of the combined 
signal. 

When noise is added in this simple case with 12 dB SNR as shown in Figure 28 
and Figure 29 there is also improved performance with a priori knowledge. 
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Figure 28 Generated Signal for A Priori Example with Noise 



los & xkkps(1 .:) (estimate = 




Figure 29 A Priori Knowledge with noise Residts 
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G. EFFECT OF Q MATRIX ON KALMAN GAINS 

The random process forcing input covariance matrix Q is a measure of the 
expected excitation of the noise. It is used in the calculation of the Kalman gains as 
shown in Chapter IV, Section A-3 on page 26 which are restated here: 

G(k) = P{ k | k - 1) • C T ■ [C ■ P(k\ k - 1) • C T + /?]"' 

P(k\k) = [/- G(k) • C] P(k\k - 1) 

P{ k + \\k) = (j>-P{k\k)-<j> T +0 



The Q matrix is a function of the state noise coefficient matrix, T , and the 
variance of the forcing function, w, as shown below. 



dt 2 




dt 



Q = T T r -var|>] 



dt 4 




2 




2 

dt 2 



• var[w] 



As there is no forcing function, w, the value of Q was initially set to zero, however 
in order to determine its impact, different values for the variance of w, var[w], were used. 
The filter results and Kalman gain values G1 and G2 are shown in Figure 30 through 
Figure 33 for different values of var[w] and Q. As the value for var[w] is increased, the 
steady state gain increases. The only noticeable effect occurs when var[w] is greater than 
100. With an increase in var[w] and Q, there is also an increase in the Kalman gains. This 
results in more responsiveness. If the gain values are too large the filter will be too 
responsive resulting in more errors. 
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Kalman Filter - Signal Options 




Figure 30 



Effect of Gain Values with varfwj=0, 



Q = 



0 0 
0 0 
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Kalman Filter - Signal Options 




Figure 31 



Effect of Gain Values with var[w]=10, 0 = 



0.0003 

0.005 



0.005 

0.1 
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Kalman Filter - Signal Options 




Figure 32 



Effect of Gain Values with var[w]=100, 0 = 



0.0025 0.05 
0.05 1 
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Kalman Filter - Signal Options 




Figure 33 



Effect of Gain Values with var[w]=150, 0 = 



0.0038 0.075 
0.075 1.5 
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Steady State Gain Values 




Figure 34 Steady State Kalman Gains with var[w]=J50 



The steady state gain values after time t=2 seconds are shown in Figure 34. For all 
other values of var[w] shown, the corresponding steady state Gain values approached 
zero. 

H. SUMMARY 

This chapter has presented the results of using the Kalman filter to estimate the 
phase in BPSK, QPSK, and Unbalanced QPSK signals to simulate the various formats of a 
GPS signal. Results have been shown for the performance with and without Gaussian 
white noise and with multipath interference. 

The designed filter works best with the BPSK signal. Gaussian white noise 
degrades performance due to the effects of the correct and incorrect phase having the 
same value at a zero crossing. It may be possible to modify the phase switching logic to 
reduce this effect by allowing a phase change to only occur at or near the expected time of 
the end of a bit period or when the error differences are greater than some minimum value. 
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The multipath interference was shown to result in errors in the estimate’s phase 
and amplitude. This error is caused by the filter tracking the phase of the combined 
multipath signal. It was shown that with a priori knowledge of multipath delay and 
amplitude, the filter is able to discriminate between the direct path and the delayed signals 
and to eliminate this error. Sources of the a priori information might include ground 
testing by placing a receiver into a self-survey mode and allowing the GPS satellites to 
pass overhead. (Lightsey) It may also be possible to use the technique suggested by 
Kumar and Lau to estimate the multipath parameters through seismic deconvolution 

For this case it was assumed that the multipath environment was not changing. In 
reality the multipath environment would be very dynamic as a result of the satellites orbit 
and changing attitude. Additionally, the attitude determination and control system only 
uses 4 out of possibly 1 1 visible satellites. The satellites selected would also effect the 
multipath environment. To overcome this may require frequent re-initialization of the 
Kalman filter. 
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VL CONCLUSION 



A. SUMMARY 

This study has focused on the carrier phase signal which is used by GPS based 
attitude determination systems. A GPS carrier signal and a multipath environment were 
described and modeled. A Kalman filter was developed to track the GPS signal and 
results were shown for Gaussian noise and in the presence of multipath. 

The developed Kalman filter was shown to work best for the L2, BPSK signal due 
to the minimum complexity involved in the phase changing logic. It was shown that the 
filter is capable of discriminating between the line of sight and delayed signals and 
performed well in the presence of multipath, given good a priori knowledge. 

More refinement is required in the phase change detection logic to eliminate errors 
caused by the error terms having the same value at a zero crossing. The Q excitation 
matrix needs considerable adjustment to handle changing phase in the presence of 
multipath. The source of the a priori information needs to be developed. 

B. OTHER METHODS OF MULTIPATH MITIGATION FOR POSSIBLE 

STUDY 

As multipath accounts for 90% of the error in an attitude determination system, 
much work is currently being done in this area. Three categories of multipath mitigation 
methods that have been used include: 1) Methods involving the geometry of the 
multipath, 2) Analyzing the consistence of the different code measurements and 3) 
Exploiting information coming from different channels. 

Of these three, the easiest to implement involves geometry. One approach is to 
focus on improving the antenna or its position. This may be accomplished by placing the 
antenna in isolated locations minimizing the multipath, something that is not always 
possible on a satellite. Another method of multipath rejection is achieved by reducing the 
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left-hand circularly polarized gain (LHCP) of the antenna without reducing the right-hand 
circularly polarized gain (RHCP). This is because the satellite signals are RHCP and the 
reflected multipath signals tend to be either linearly polarized or even LHCP, depending 
on the reflecting surface. (Van Dierendonck) An example of an antenna designed 
specifically to reduce multipath effects is the choke ring antenna. This antenna is 
comprised of vertical-aligned concentric rings centered about the antenna element that are 
connected to the ground plane. The vertical rings shape the antenna pattern such that the 
multipath signals incident on the antenna at the horizon and negative elevation angles are 
attenuated (Kaplan). 

Other methods of dealing with the multipath problem involve modification to the 
GPS receiver for code measurements. One approach involves reducing the early-late 
delay spacing among the correlators in the GPS receiver code lock loop. This technique 
provides limited results if the early-late spacing is smaller than the initial delay error due to 
multipath. (Kumar) 

An example of the third category is the method presented by Kumar and Lau 
involving the use of optimal deconvolution to estimate the impulse response of the 
effective multipath channel and obtain an inverse filter which equalizes the multipath 
channel (Kumar) Another technique involves adaptively estimating the spectral 
parameters (frequency, amplitude, phase offset) of multipath in the associated SNR and 
then constructing a profile of the multipath error in the carrier phase. A multipath 
correction is then made by subtracting the profile from the actual phase measurement. 
(Comp) 
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APPENDIX A 



MATLAB PROGRAMS 



A. DOKAL.M 



% File Name: do kal.m 
% Last Updated: 1 8 Apr 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Script file to set up and run kalman filter simulations 

clear 

uivars 

% End of File do_kal.m 



B. Gl.M 



% File Name: gl.m 
% Last Updated: 18 Oct 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function generates the Gl(t) signal 
% for GPS C/A Code generator 

% 

% 



function [new R, g] = gl(R, xl) 



if xl==l 

R = ones(size(R)); % Reset to all ones at XI epoch 
end; 

new_R=gen_poly([l,3, 10], R); 
g=R(10); 



% End of File gl.m 
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c. 



G2.M 



% File Name: g2.m 
% Last Updated: 1 8 Oct 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function generates the G2(t) signal 
% for GPS C/A Code generator 
% 

% 



function [newR, g] = g2(R, xl, sv) 



if xl==l 

R = ones(size(R)); % Reset to all ones at XI epoch 
end; 

new_R = gen_poly([l 2 3 6 8 9 10], R); 

%new_R(l) = mod2(R(2)+R(3)+R(6)+R(8)+R(9)+R( 10)); 
%new_R(2:10) = R(l:9); 



% phase select logic, see p. 92 KAPLAN 
% This changes based on satellite choosen. 
sv_tap =[26 
37 
48 
59 

1 9 

2 10 
1 8 
29 
3 10 
23 
34 
56 
67 
78 

8 9; 9 10;1 4;2 5;3 6; 4 7; 5 8; 6 9; 1 3; 4 6; 5 7; 6 8; 7 9; 
8 10; 1 6; 2 7; 3 8; 4 9; 5 10; 4 10; 1 7; 2 8; 4 10]; 



g=xor( R(sv_tap(sv,l) ), R(sv_tap(sv,2)) ); 



% End of File g2.m 
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D. 



GEN POLY.M 



% File Name: gen_poly.m 
% Last Updated: 18 Oct 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function generates the linear code generators 
% the polynomial is of the form 1 + sum(x(i)) where 
% sum(x(i)) means that the output of the i'th cell of the register 
% is used as the input to the modulo-2 adder (exclusive or) 

% and the 1 means that the output of the adder is fed to the first cell 
% 

% so 1 + xl + x3 + xlO could be generated with poly =[113 10] and R 
% is the register to perform the operation on 

function r = gen_poly(poly, R) 

result = mod2( sum(R(poly(2:length(poly) ) ))); % xor registers according to 

% polynomial 

r(2:length(R)) = R(1 :(length(R)-l)); % shift registers right 

r(poly(l))=result; % put result in register according 

% to polynomial 



% End of File gen_poly.m 



E. INIT.M 



% File Name: init.m 
% Last Updated: 1 1 Nov 96 
% Written By: T. H. Newman 
% Operating Environment: Sun Openwin & Win 95 
% Matlab Version 4.2c. 1 
% 

% Description: This script file sets the parameters for 
% the GPS signal simulation 

global w; 

% my _glb; 
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fo = 10.22999999 543e6; 

fl = 120*fo; % frequency of L2 signal 

f2 = 154*fo; % frequency of LI signal 



To = .001; 

codelength =10; 
fsl.s2,s31=si g g enfcode length); 
sigl =[( 1 :length(sl))',sT]; 
sig2= [( 1 : Iength(s2))',s2'] ; 
sig3=[(l :length(s3))',s3']; 

kalinit; 

% End of File init.m 



F. KALINIT.M 



% File Name: kal_init.m 
% Last Updated: 27 Feb 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Initialization file for Kalman Filter 

% 



global Pkkml; 
global R; 
global Q; 
global I; 
global Phi; 
global xkkml; 
global Tb; 
global Noise _gain; 
global dt; 

global w; 

% my_glb; 

A=[0 1; -w A 2 0]; % xl = cos(wt) 

% x2 = xldot = -w*sin(wt) 
% x2_dot = -w A 2*cos(wt) 

% x = [xl x2]' 
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B=[0 w A 2]’; 



% x_dot = A*[xl_dot, x2_dot] 
% => A= [0 l;-w A 2 0]’ 

% x_dot = Ax + Bu 
% u is control input 



C=[l 0]; % will use to get xl estimate 

D=0; 



I=eye(2); % Identity Matrix 

[Phi,Del]=c2d(A,B,dt); % converts the continuous-time system: 

% x_dot = Ax + Bu 
% to the discrete-time state-space 
% sy stem: 

% x[n+l] = Phi * x[n] + Del * u[n] 
% assuming a zero-order hold on the 
% inputs and sample time dt 



Pkkml=le6*I; 



R=l; 



%Q=0*I; 



% Initialize P(k|k-1) =M 

% where M = covariance of initial state 
% want this number large to give 
% weight to 1st measurment 
% var[e(0|0)] = var[v(0)] = R 
% v is reference input 
% (noise in measurement) 

% Q is measure of exitation by forcing function 



var_w=150; 

Q = [ (dt A 4)/4 (dt A 3)/2; (dt A 3)/2 dt A 2 ]*var_w; 



xkk = [0 0; 0 0]; % initialize x(k|k) to all zeros 

% 



xkkps= [0 0; 0 0]; 



x(k|k) is state of plant at time k given k observ ations 



xkkml = [0,0; 0,0]; 



% initialize x(kjk-l) => x(0|-l) = x_o = 0 



% 

% clear all the variables used for saving plotting information 

clear myrand; 

clear x; 

clear sm; 

clear x_delay; 

clear z; 

clear y; 

clear xkk__ps; 

clear t; 
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% End of File kal init.m 



G. KALMAN.M 



% File Name: kalman.m 
% Last Updated: 28 Feb 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Kalman filter for tracking phase of a 
% QPSK signal 
% 

% Based on Professor Titus’ algorithim which follows Prof. Kirk's 
% NPS lecture notes " Optimal Estimation: An Introduction to the 
% Theory and Applications" 1975 All equation numbers are referenced 
% to the lecture notes. 

function y=kalman(z) 

global Pkkml ; % Values initialized in kal_init 

global R; 

global Q; 

global I; 

global Phi; 

global xkkml; 

global w; 



C=[1,0]; 
Cps = [0,1]; 



G=Pkkml*C'*inv(C*Pkkml*C' +R); 



% Calculate Gain at time k 



% G(k) = P(k|k- 



l)*C'*inv(C*P(k|k-l)*C' + R) 



% cqn 4.3-102a 



Pkk=(I-G*C)*Pkkm 1 ; % Calculate P(k|k) 



% P(k|k) = [I -G(k)*C(k)]*P(k|k-l) eqn 4.3-103a 



Pkkm 1 =Phi*Pkk*Phi' + Q; % P(k|k- 1 ) 



eO = z-C*xkkml; 



xO=xkkml+G*eO; 



el = z-Cps*xkkml/w r ; xl=Cps*(xkkml+G*el)Av; 



e2 = z+C*xkkml; 
e3 = z+Cps*xkkml/w; 



x2=-(xkkm 1 +G*e2); 
x3= -Cps*(xkkml+G*e3)/w; 



if abs(eO)<=abs(el); 
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e=eO; % Smallest error is with 

xkk = xkkml +G*(e); % old phase 

xkkps = xkk; 
else; 

e=el ; % Smallest error is +90 

xkk = xkkml +G*(e); % is 
xkkps = Cps*xkk/w; 
end 

if abs(e2) <= abs(e) 
e=e2; 

xkk - xkkml + G*(e); 

xkkps = -xkk; % Smallest error is -90 

end; 

if abs(e3) <= abs(e) 
e=e3; 

xkk = xkkml + G*e; % Smallest error is 180 

xkkps = -Cps*xkkAv; 
end; 



xkkml =Phi*xkk; % eqn 4.3-89a 

y 1= C*xkk; % report the non shifted version 

y=[vl(U); xkkps(l); x0(l); xl(l); x2(l); x3(l)]; 



% End of File kalman.m 



H. KALMAN U.M 



% File Name: kalman_u.m 
% Last Updated: 26 May 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Kalman filter for tracking phase of a 
% unbalanced QPSK signal 

% 

% Based on Professor Titus’ algorithim which follows Prof. Kirk's 
% NPS lecture notes ” Optimal Estimation: An Introduction to the 
% Theory and Applications" 1975 All equation numbers are referenced 
% to the lecture notes. 
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function y=kalman(z) 



global Pkkm 1 ; % Values initialized in kal_init 

global R; 

global Q; 

global I; 

global Phi; 

global xkknil; 

global w; 

global Al; 

global A2; 

011 , 0 ]; 

Cps = [0,1]; 

G=Pkkm 1 *C'*inv(C*Pkkml *C' +R); 
l)*C'*inv(C*P(k|k-l)*C' + R) 



Pkk=(I-G*C)*Pkkm 1 ; % Calculate P(k|k) 

% 

Pkkm 1 =Phi*Pkk*Phi' + Q; % P(k|k-1) 



% Calculate Gain at time k 

% G(k) = P(k|k- 

% eqn 4.3-102a 

P(k|k) = [I -G(k)*C(k)]*P(k|k-l) eqn 4.3-103a 



F=( 1/A2); % Factor for imbalanced quadature phase term 

cos_\vt=C*xkkml ; 
sin_wt=-Cps*xkkm 1/w; 



eO = z-(cos_wt + sqrt(2)*sin_wt); 
el = z-(cos_wt - sqrt(2)*sin_wt); 
e2 = z+(cos_wt - sqrt(2)*sin_wt); 
e3 = z+(cos_wt + sqrt(2)*sin_wt); 

if abs(eO)<=abs(el); 

e=eO; % Smallest error is with 

xkk = xkkml +G*(e); % old phase 

xkkps = (cos_wt+sqrt(2)*sin_wt); 
else; 

e=el; % Smallest error is +90 

xkk = xkkml +G*(e); % is 
xkkps = ( cos_wt - sqrt(2)*sin_wt ); 
end 



if abs(e2) <= abs(e) 
e=e2; 

xkk = xkkml + G*(e); 

xkkps = (cos_wt - sqrt(2)*sin_wt); % Smallest error is -90 

end; 

if abs(e3) <= abs(e) 
e=e3; 
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xkk = xkkml + G*e; % Smallest error is 180 

xkkps = (cos_vvt + sqrt(2)*sin_wt); 



end; 



xkkml =Phi*xkk; 



% eqn 4.3-89a 



yl= C*xkk; 



% report the non shifted version 



y=[yl(l,l); xkkps(l)]; 



% End of File kalman_u.m 



I. KALMAN2.M 



% File Name: kalman2.m 
% Last Updated: 3 May 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Kalman filter for tracking phase of a 
% BPSK signal. 

% 

% Based on Professor Titus' algorithm which follows Prof. Kirk’s 
% NPS lecture notes " Optimal Estimation: An Introduction to the 
% Theory and Applications" 1975 All equation numbers are referenced 
% to the lecture notes. 

function y=kalman(z) 

global Pkkml ; % Values initialized in kal_init 

global R; 

global Q; 

global I; 

global Phi; 

global xkkml; 

global w; 



C=[1,0]; 
Cps = [0,1]; 



G=Pkkml *C'*inv(C*Pkkml *C' +R); 



% Calculate Gain at time k 



% G(k) = P(k|k- 



1 )*C'*inv(C*P(k|k- 1)*C* + R) 



% eqn4.3-102a 
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Pkk=(I-G*C)*Pkkml; % Calculate P(k|k) 

% P(k|k) = [I-G(k)*C(k)]*P(k|k-l) eqn 4.3-103a 

Pkkml=Phi*Pkk*Phi' + Q; % P(k|k-1) 



eO = z-C*xkkml; xO=xkkml+G*(eO); 
e2 = z+C*xkkml; x2=-(xkkml+G*e2); 



if abs(eO)<= abs(e2); 

e=eO; % Smallest error is with 

xkk = xkkml +G*(e); % old phase 

xkkps = xkk; 
else 
e=e2; 

xkk= xkkml + G*(e); 

xkkps = -xkk; % Smallest error is -90 

end; 



xkkml=Phi*xkk; %eqn4.3-89a 

yl= C*xkk; % report the non shifted version 

y=[yl(l,l); xkkps(l)]; 



% End of File kalman2.m 



J. MOD2.M 



% File Name: mod2.m 
% Last Updated: 18 Oct 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Returns the modulo 2 value of input 
% if input is multiple of 2 returns 0 
% else returns 1 

% 

% 



function x = mod2(y) 



if (y/2)* 10 == round(y/2)* 10 
x = 0; 
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else 

x=l; 

end 



% End of File mod2.m 



K. MPATH.M 



% File Name: mpath.m 
% Last Updated: 28 Apr 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function simulates the multipath environment 

function sm = mpath(s) 

global delay 
global alpha 
global tau 

len=length(delay); 

if ( len > 1 ) 

delay (2: len) = delay(l :len-l); % Shift delay vector 
delay(l) = s; 

sm = sum(delay(tau).*alpha); 
else 

delay(l)=s; 

sm=s; 

end; 



% End of File mpath.m 



L. MPATHAP.M 



% File Name: mpathap.m 
% Last Updated: 4 Jun 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 
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% Description: Kalman filter for multipath with a prion 
% information. 

clear; % clear old variables 

w=2; % = 2*pi*f 

Tf= 15; % finish time ~ sec 

dt=.l; 

noise_fac=0; % change scale on noise 

%multipath a priori information 

alpha = .5; % delayed signal amplitude 

tau = 10; % time of delay 

Amp = 1 ; % amplitude of line of sight signal 



A=[0 1 0 0; -w A 2 0 0 0;0 0 0 1;0 0 -w A 2 0]; 
B=[0 w’ A 2 0 0]'; 

C=[l 0 1 0], 

[Phi,Del]=c2d(A,B,dt); 

I=eye(4); 

Pkkml=le6*I; 

R=l; 

Q=0.1*I; 
kmax=Tf/dt +1 ; 



switch 1 = 3; % time of phase changes 

switch2 = 9; 

time=0; 

for (i=l :kmax) % generate line of sight signal 

t(i)=time; 
if t(i) <= switch 1 

los(i) = cos(w*t(i)); % line of sight direct path signal 
elseif t(i) >= switch 1 & t(i) <= switch2 
los(i) = sin(w*t(i)); 
else 

los(i) = cos(w*t(i)); 
end; 

time=time+dt; 
end; 
tl=t; 
time = t; 
clear t; 

for i= 1 :length(los)-tau % generate delayed signal 

x(i) = los(i+tau); % Line of sight signal x(t) 

delay(i) = los(i)*alpha; % delayed signal alpha*x(t-tau) 
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z(i) = x(i) + delay(i); % z = x(t) + alpha*x(t-tau) 
t(i) = time(i+tau);; 
end; 

z = z+noise_fac*randn(size(z)); % add in noise 

\vt = w*time(tau+l); 
wtmtau = w*time(l); 
xkkml=[ cos(wt) 

-w*sin(wt) 
alpha*cos( wtmtau) 

-w*alpha*sin( wtmtau)] ; 



for (i=l:length(z)) 

G=Pkkm 1 * C’*inv(C*Pkkm 1 *C' +R); 
Pkk=(I-G*C)*Pkkml; 
Pkkml=Phi*Pkk*Phi' + Q; 
eO=z(l,i)-[l 0 1 0]*xkkml(:,i); 

el=z(l,i)-[0 -1/w 1 0]*xkkml(:,i); 

e2=z(l,i)-[0 -1/w 0 -l/w]*xkkml(:,i); 
e3=z(l,i)-[l 0 0 -l/w]*xkkml(:,i); 
%e_vect(:,i) = [ eO, el, e2, e3 ]'; 



if abs(eO)<=abs(el) & abs(eO) <= abs(e2) % abs(e0)<=abs(e3); % t(i)<=sl; 

%if (t(i) <= switch 1) 
e=eO; 

xkk(:,i) = xkkml(:,i) +G*(e); 
xklqK(l,i) = xkk(l,i); 
xkkps(3,i) = xkk(3,i); 

elseif abs(el) <= abs(eO) & abs(el) <= abs(e2) & abs(el)<=abs(e3) % si < t < sl+phi 
%elseif (switchl <= t(i) & t(i) <= switchl+tau*dt) 
e=e 1 ; 

xkk( : ,i)=xkkm 1 (: ,i)+G*(e); 
xkkps(l,i) = -xkk(2,i)/w; 
xkkps(3,i) = xkk(3,i); 

elseif abs(e2)<=abs(e0) & abs(e2)<=abs(el) & abs(e2)<=abs(e3) % sl+tau < t < s2 
%elseif (switch l+tau*dt <= t(i) & t(i) <= switch2) 
e=e2; 

xkk(:,i)=xkkml(:,i)+G*e; 
xkkps(l,i) = -xkk(2,i)/w; 
xkkps(3,i) = -xkk(4,i)/w; 

elseif abs(e3)<=abs(e0) & abs(e3)<=abs(el) & abs(e3)<=abs(e2) 

%elseif (switch2 <= t(i) & t(i) <= switch2+tau*dt) 

% s2 < t < s2+phi 

e=e3; 

xkk(: ,i)=xkkm 1 (:,i)+G*e; 
xkkps( 1 ,i)=xkk( 1 ,i); 
xkkps(3,i)=-xkk(4,i)/w; 
else 

e=eO; 
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xkk(:,i) = xkkml(:,i) +G*(e); 
xkkps(l,i) = xkk(l,i); 
xkkps(3,i) = xkk(3,i); 

end; 

xkkml (:,i+l)=Phi*xkk(:,i); 

end; 



%%%%%%%%%%%%%%%%%%%%% 0 /o 0 /o 0 /o 0 /o%% 0 /o% 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o 0 /o% 0 /o 0 /o 0 /o% 0 /o%%%% 

%%%%%%% 



figure; 

mpath=delay; 

los=x; 

figure; 

subplot(3,l,l); 

plot(t,los,'r'); 

ylabel('los'); 

subplot(3,l,2); 

plot(t,mpath,'b’); 

ylabel('mpath'); 

subplot(3,l,3); 

plot(t,z,'g'); 

ylabel('z'); 

figure; 

subplot(2,l,l) 

plot(t,los,'r',t,xkkps(l,:),tK'); 
title('los & xkkps(l,:) (estimate = 

subplot(2,l,2); 
plot(t,mpath,'r',t,xkkps(3 , 
title('mpath & xkkps(3,:)'); 



% End of File mpath_ap.m 
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M. 



MY GLB.M 



% File Name: my_glb.m 
% Last Updated: 16 Feb 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Global declaration for gps simulation and 
% kalman filter. 



global C; % Values initialized in kal_init.m 

global Pkkml; 
global R; 
global Q; 
global I; 
global Phi; 
global xkkml; 

global w; 

fl = 1; 

To = 1; 
w=2*pi*fl; 
phase_init = pi/6; 
num_save = 100000; 



% NOTE: These are the same as in init.m 
% and are placed here for use in non-gps model. 

% Set frequency for simulation 
% Phase offset for sine wave in BPSK modulator 
% Number of points to save to workspace for plotting later 



% End of File my_glb.m 



N. SHIFTREG.M 



% File Name: shift_reg.m 
% Last Updated: 21 Oct 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function generates the shift register 
% linear code generators where R = the input values 
% size = register size. The register is shifted right 
% and returned in a new vector 

function r = shift_reg(R, size) 

r = ones(Lsize) 



% End of File shift_reg.m 
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o. 



SIG GEN.M 



% File Name: sig_gen.m 
% Last Updated: 11 Nov 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Returns either 1 P(Y) code xor with data 
% 2 P(Y) code 

% or 3 C/A code xor with data 

% the number is the sig num which determines which signal 
% is returned 
% 



function [sl,s2,s3] = sig_£en(code_length) 
sv=l; 

xlepoch = 0; 
xlA_epoch = 0; 
xlB_epoch = 0; 
x2A_epoch = 0; 
x2B_epoch = 0; 

%code_length = 10; 

r_gl = ones( size(l:10) ) ; % initialize register values 

r_g2 = ones( size(l : 10) ); 

r_xla =[0 0010010010 0]; 

rxlb =[001010101010]; 

r_x2a = [10100100100 1]; 

r_x2b =[001010101010]; 

shiftregister = zeros(size(l:37)); 

t = 0; 

data_clock = 20; 



for i = 1 .codelength 
if data_clock == 20 

data =1; % set to all ones to ignore 

data_clock = 0; 

else 



end; 



dataclock = dataclock + 1 ; 



forj =1:10 

[r_gl, g_l] = gl(r_gl, xl_epoch); 



% gen C/A code at 1 .023 MHz 
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[r_g2, g 2] = g2(r_g2, xl_epoch, sv); 

g = xor(g_l, g_2); 

for k=l:10 % gen P/A code at 10.23 MHz 

t=t+l; % or 10 * C/A code generation rate 

[r_xla, A] = xla(r_xla, xlAepoch); 

[r xlb, B] = xlb(r_xlb, xlB_epoch); 

[r_x2a, C] = x2a(r_x2a, x2A_epoch); 

[r_x2b, D] = x2b(r_x2b, x2B_epoch); 

xl=xor(A.B); 

x2=xor(C,D); 

shift_register(2:37) = shift_register(l:36); 
shift_register(l) = x2; 
p=xor(xl, shiftregister(sv)); 

sl(t) = xor(p,data); 
s2(t) = p; 

s3(t) = xor(g. data); 
end; % for k; 

end; %forj 

end; % for i 



% End of File sig gen.m 



P. TEST_CA.M 



% File Name: test_CA.m 
% Last Updated: 31 Mar 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Script file to test C/A code generation. Output 
% should match First 10 C/A-Chips (Octal) of GPS-ICD-2000 table 

for sv=l:37 

r_gl = round( rand( size(l : 10) ) ); % initialize to some random state 
r_g2 = round( rand( size(l : 10) ) ); 

xl = 1; 

i=i; 

[ r _gk g_l(i)] = gl(r _gl, xl); % Initialize at epoch 

[r_g2, g_2(i)] = g2(r_g2, xl, sv); 
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% Simulate 1st 10 chips 



xl=0; 
for i=2:10 

[r_gl, g_l(i)] = gl(r_gl, xl); 
[r_g2, g_2(i)] = g2(r_g2, xl, sv); 
end; 



g = xor(g_l, g_2); 



% 



digl = g(l); 

dig2 = g(2)*4 + g(3)*2 + g(4); % convert to octal number 

dig3 = g(5)*4 + g(6)*2 + g(7); 

dig4 = g(8)*4 + g(9)*2 + g(10); 

result(sv) = digl*10 A 3 + dig2*10 A 2 + dig3*10 + dig4; 

end 

result % display results 

% End of File test CAm 



Q. TEST_P.M 



% File Name: test_P.m 
% Last Updated: 1 April 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Script file to test P code generation. Output 
% should match First 12 P-Chips (Octal) of GPS-ICD-2000 table 
% Note that only the first 13 sv’s are calculated. SV 9 - 37 are 
% identical 

for sv =1:13 

r_xla = round( rand( size(l:12) ) ); % initialize to some random values 
r_xlb = round( rand( size(l:12) ) ); 
r_x2a = round( rand( size(l:12) ) ); 
r_x2b = round( rand( size(l : 12) ) ); 

epoch = 1; 

i=l; 

[r_xla, A(i)] = xla(r_xla, epoch); 

[r xlb, B(i)] = xlb(r_xlb, epoch); 

[r_x2a, C(i)] = x2a(r_x2a, epoch); 

[r_x2b, D(i)] = x2b(r_x2b, epoch); 



epoch=0; 
for i=2: 12 

[r xla, A(i)] = xla(r_xla, epoch); 
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[r_xlb, B(i)J = xlbtrxlb, epoch); 

[r_x2a, C(i>] = x2a(r_x2a, epoch); 

[r_x2b, D(i)] = x2b(r_x2b, epoch); 
end; 

xl=xor(A,B); 

x2=xor(C,D); 

delay=sv; 

xl=[xl, zeros(l, delay)]; 
x2=[ones(l, delay), x2 j; 
p=xor(xl, x2); 

digl = p(l)*4 + p(2)*2 + p(3); 
dig2 = p(4)*4 + p(5)*2 + p(6); 
dig3 = p(7)*4 + p(8)*2 + p(9); 
dig4 = p(10)*4 + p(ll)*2 + p(12); 
result(sv) = digl*10 A 3 + dig2*10 A 2 + dig3*10 + dig4; 
end; % for sv 
result 

% End of File test P.m 



R. UIINIT.M 



% File Name: uiinitm 
% Last Updated: 4 May 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: Called by ui_vars.m to initialize parameters 
% after the Init button is pushed. 

% 



delay=zeros(size(l:max(tau) )); 
if SNR >0 

Noise_gain = sqrt((Amp A 2)/( 10 A (SNR/10) )); 
else 

Noise_gain = 1; 
end; 

if SNR >=100 
Noisegain = 0; 
end; 

w=2*pi*f; 

disp(sprintf(' f=%2.2f Initial Phase = %2.2f , f, p degree)) 
disp(sprintf(' Amp=%2.2f, Amp)) 
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disjKsprintfC dt=%2.2f, stop = %2.2f, Tb=%2.2f , dt, stop, Tb» 
disp(sprintf(' SNR=%2.2f dB, Noise_gain=%2.2f ,SNR, Noise_gain)) 
disp(’tau='), disp(tau) 
disp('alpha='),disp(alpha) 



% End of File ui init.m 



S. UIVARS.M 



% File Name: ui_vars.m 
% Last Updated: 28 Apr 97 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: An Editable Text windows to control 
% initialization values for kal init 
% Calls: my_global.m to initialize values 
% kal_init.m to initialize kalman filter variables 

% init.m to generate GPS code sample 

global MainFig 

% below are the initial values that will be in the text window 

global C; % Values initialized in kal_init.m 

global Pkkml; 

global R; 

global Q; 

global I; 

global Phi; 

global xkkml; 

global Tb; 

global Noise_gain; 

global dt; 

global Amp; % Amplitude for BPSK 

global tau; 

global alpha; 

global delay; 

global w; 

global clock; 

my_glb; 

% This initilizes w, phase init, num_save, Tb, Noise gain, Al, A2, dt 
% Retreive initial values from the strings 
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stan = 0; 
stop = 10: 
f = w/(2*pi); 



%stan time 



%stop time 
%signal frequency 



p_degree = phase_init* 180/pi; %phase in degrees 
SNR = 100; 

f_string = sprintfC %2.2f ,f);; 

p_degree_string = sprintf(' %2.2f ,p_degree); 

Amp_string = sprintfC %2.2f ,Amp); 

dt_string = sprintfC %2.2f ,dt);; 

start_string = sprintfC %2.2f, start); 

stop_string = sprintfC %2.2f ,stop); 

Tb_string = sprintfC %2.2f .Tb); 

SNR string = sprintfC %2.2f, SNR); 

tau_string = ' [1]'; 

alpha_string = ' [1]’; 

tau=eval(tau_string) ; 
alpha = eval(alpha_string); 

%define the Main Figure Window. 

numentries = 9; 

sep=45; 

ymax = sep*(num_entries)+ 1 0; 

x_pos=25; y_pos=35; dx=125; dy = ymax; 

MainFig = figure ('Position', [x_pos y_pos dx dy],,.. 

'Color', 'white', 'NumberTitle', 'off, 'Name',... 

'Simulation Variables',... 

'MenuBar','none'); 

% % 

% Define UIMENU for pull down menu options 

e=uimenu('laber,'Init', 'callback', . . . 

1cal_init,ui_init; '); 

f=uimenu(’laber, 'Plot'); 

uimenu(f, 'label', 'Phase Lock Loop', 'callback', 'plot_k'); 
uimenu(f, 'label', 'Delay', 'callback', 'plotd'); 
uimenu(f, 'label', 'QPSK','callback','plot_q2'); 
uimenu(f, 'label', 'BPSK', 'callback', 'plotbpsk'); 
uimenu(f,'laber,'Error','callback','plot_err'); 
uimenu(f, 'label', 'Variables', 'callback', ’plot_k2'); 

g=uimenu('label', 'Sim'); 
gl=uimenu(g, ’label','BPSK’); 
uimenu(gl, 'label', 'BPSK Simulation', 'callback', 'kal sim'); 
uimenu(gl, 'label', 'PLL Simulation','callback’,'pH’); 
g2=uimenu(g, 'label', 'QPSK'); 
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uimenu(g2, 'label', 'QPSK Simulation’, 'callback','kalq_sim'); 
uimenu(g2, 'label', 'PLL Simulation', 'callback', 'pll_q'); 
g3=uimenu(g, 'label', 'Unbalanced QPSK'); 
uimenu(g3, 'label', 'Unbalanced QPSK Simulation', 'callback', 'kalusim'); 
uimenu(g3, 'label', PLL Simulation', 'callback', 'pll u'); 
g4=uimenu(g, 'label', 'Multipath', 'callback','mpathsim'); 



n=l; %uimenu number 

ymax = dy; 

xpos=10; 

dx = 100; 

dy-20; 

space = 22; 

label_fg_color = 'Black'; 
label_bg_color = 'White'; 
box_fg_color = 'Black'; 
box_bg_color = 'Yellow'; 

f_label = uicontrol(MainFig, 'Style', 'text', 'Pos',[xpos ymax-n*sepH-space dx dy],... 
'String', 'Signal freq (Hz)',ToregroundColor',label_fg_color,... 
*BackgroundColor',label_bg_color); 



f_box = uicontrol(MainFig, 'Style', 'edit', 'String', f_string,... 

'Pos',[xpos ymax-n*sep dx dy],... 

Toregroundcolor',box_fg_color,.. . 

'BackgroundColor',box_bg_color, . . . 

'Callback','f = eval( get(f_box,"String") ); w=2*pi*f;'); 

n=n+l; 

phasejabel = uicontrol(MainFig, 'Style', 'text', ’Pos',[xpos ymax-n*sefri-space dx dy],... 
'String', 'Phase Offset °',’ForegroundColor',label_fg_color,... 

'BackgroundColor' , labelbgcol or) ; 

phasebox = uicontrol(MainFig, 'Style', 'edit', 'String', pdegreestring,... 

'Pos’,[xpos ymax-n*sep dx dy],... 

Toregroundcolor',box_fg_color, . . . 

'BackgroundColor', box_bg_color, . . . 

'Callback', 'p degree = eval( get(phase_box, "String”) );phase_init=p_degree*pi/180'); 
n=n+l; 

amp label = uicontrol(MainFig, 'Style', 'text', 'Pos',[xpos \max-n * sep+space dx dy],... 
'String', 'Amplitude', 'ForegroundColor',label_fg_color, . . . 

'BackgroundColor', labelbgcolor); 

amp box = uicontrol(MainFig,'Style','edit', 'String', Amp string,... 

'Pos',[xpos ymax-n*sep dx dy],... 
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'Foregroundcolor' ,box_fg_color, . . . 
'BackgroundColor',box_bg_color, . . . 

'Callback', 'A1 = eval( get(amp_box,"String") );'); 



n=n+l; 

samp_label = uicontrol(MainFig, 'Style','text', 'Pos',[xpos ymax-n*sep+space dx dy], 
'String', 'Ts (sec)',ToregroundColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 

samp_box = uicontrol(MainFig, 'Style', 'edit', 'String', dt_string,... 

■Pos',[xpos ymax-n*sep dx dy'],... 

'Foregroundcolor',box_fg_color,. . . 

'BackgroundColor',box_bg_color, . . . 

'Callback', 'dt = eval( get(samp_box,"String") );'); 

n=n+l; 

stop_label = uicontrol(MainFig, 'Style', 'text', 'Pos',[xpos ymax-n*sep+space dx dy],.. 
'String', 'Stop time (sec)','ForegroundColor',label_fg_color,... 
’BackgroundColor',label_bg_color); 

stop_box = uicontrol(MainFig,’Style', 'edit', 'String', stop_string,... 

'Pos',[xpos ymax-n*sep dx dy],... 
r Foregroundcolor’,box_fg_color, . . . 

'BackgroundColor',box_bg_color, . . . 

’Callback', 'stop = eval ( get(stop_box, "String") 

n=n+l; 

Tb_label = uicontrol(MainFig, 'Sty'le','text', *Pos',[xpos ymax-n*sep+space dx dy],... 
'String', 'Tb (sec)','ForegroundColor',label_fg_color,... 
’BackgroimdColor',label_bg_color); 

Tb_box = uicontrol(MainFig,'Style', 'edit', 'String', Tb string,... 

'Pos',[xpos ymax-n*sep dx dy],... 

'Foregroundcolor', box_fg_color, . . . 

'BackgroimdColor',box_bg_color, . . . 

'Callback', 'Tb = eval ( get(Tb_box, "String") );'); 



n=n+l; 

SNR_label = uicontrol(MainFig, 'Style', 'text', 'Pos',[xpos yTnax-n*seFH-space dx dy],. 
'String', 'SNR dB','ForegroundColor',label_fg_color,... 
'BackgroundColor',label_bg_color); 

SNR_box = uicontrol(MainFig,'Style','edit','String', SNR_string,... 
'Pos',[xposymax-n*sep dx dy],... 

'Foregroundcolor',box_fg_color,. . . 

'BackgroundColor',box_bg_color, . . . 

'Callback','SNR = eval( get(SNR_box,"String") ); '); 
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n=n+l; 



tau_label = uicontrol(MainFig, 'Style', 'text', ’Pos',[xpos ymax-n*sep+space dx dy],... 
'String', 'tau','ForegroundColor',label_fg_color,... 
’BackgroundColor',label_bg_color); 

tau_box = uicontrol(MainFig, 'Style', 'edit','String', tau_string,... 

'Pos',[xpos ymax-n*sep dx dy],... 

'Foregroimdcolor',box_fg_color,. . . 

'BackgroundColor',box_bg_color, . . . 

'Callback','tau = eval( get(tau_box, "String") ); '); 

n=n+l; 



alphajabel = uicontrol(MainFig, ’Style', 'text', 'Pos',[xpos ymax-n*sep+space dx dy],... 
'String', 'alpha','ForegroundColor',label_fg_color,... 
’BackgroimdColor',label_bg_color); 

alpha_box = uicontrol(MainFig, 'Style', 'edit', 'String', alpha_string,... 

'Pos',[xpos ymax-n*sep dx dy],... 

'Foregroundcolor', box_fg_color,. . . 

'BackgroundColor' ,box_bg_color , . . . 

'Callback', 'alpha = eval( get(alpha_box,"String") ); '); 

n=n+l; 

% % 

% Define the UICONTROL button for callbacks. 

% bl =uicontrol('style', 'pushbutton', 'string', 'Initialize', . . . 

% 'position', [30 ymax-n*sep 65 25], 'Callback',... 

% 'kal init; delay=zeros(size(l:max(tau) ));'); 



% End of File ui vars.m 



T. X1A.M 



% File Name: xla.m 
% Last Updated: 1 8 Oct 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 
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% Description: This function generates the xla(t) signal 
% for GPS P Code generator 

function [new_R, A] = xla(R, xla_epoch) 

if xla_epoch==l 

R = [0 0010010010 0]; % Reset at X1A epoch 
end; 

new_R = gen_poly([l,6,8,ll,12], R); 

A=R(12); 

%EndofFilexla.m 



U. X1B.M 



% File Name: xlb.m 
% Last Updated. 18 Oct 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function generates the xlb(t) signal 
% for GPS P Code generator 

function [new_R, B] = xlb(R, xlb_epoch) 

if xlb_epoch=l 

R= [00 1 0 1 0 1 0 1 0 1 0]; % Reset at X1B epoch 
end; 

new_R = gen_poly([l 1 2 5 8 9 10 1 1 12], R); 

B=R(12); 

% End of File xlb.m 



V. X2A.M 



% File Name: x2a.m 
% Last Updated: 18 Oct 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function generates the x2a(t) signal 
% for GPS P/A Code generator 
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function [new R, C] = x2a(R, x2a_epoch) 



if x2a_epoch==l 

R = [1 0 1 0 0 1 0 0 1 0 0 1]; % Reset at X2a epoch 
end; 

new_R = gen_poly([l 1 3 4 5 7 8 9 10 11 12], R); 
C=R(12); 

% End of File x2a.m 



W. X2B.M 



% File Name: x2b.m 
% Last Updated: 18 Oct 96 
% Written By: T. H. Newman 
% Operating Environment: Windows 95 
% Matlab Version 4.2c. 1 
% 

% Description: This function generates the x2b(t) signal 
% for GPS P Code generator 

function [newJR, D] = x2b(R, x2b_epoch) 



if x2b_epoch==l 

R = [0 0 1 0 1 0 1 0 1 0 1 0]; % Reset at X2B epoch 
end; 

newJR = gen_poly([l 2 3 4 8 9 12], R); 

D=R(12); 

% End of File x2b.m 
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APPENDIX B SIM U LINK DIAGRAMS 



A. BPSK.M 



Inport 



■*JV 



Zero Order Hold 



Relay 



& 



A1*sin(wt+phi) 



Mux 



Product 



Mux 



L> 1 
Outport 



Figure 35 BPSK.M 



-> [ BPSK 
To Workspace 
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B. BPSK P90.M 



Inport 






Zero Order Hold 



Relay 



A*sin(wt + 90) 
=A*cos(wt) 



Product 



>[ Mux 

* 



BPSK_p90 



To Workspace 



Mux 



1 



Outport 



Figure 36 BPSKP90.M 
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C. BPSKU.M 




Inport 






Zero Order Hold 



sqrt(2)*A1* 

sin(wt+phi) 



Product 



-> Mux \ - k : BPSK ] 
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Mux 
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Outport 



Figure 37 BPSK U.M 



D. KALSIM.M 
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Figure 38 KAL SIMM 
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E. KALQSIM.M 



0 M ~ 

Clock time 
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Figure 39 KALQ SIM.M 



F. KALUSIM.M 

© HZ LJ 

Clock time 
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Figure 40 KALUJSIM.M 
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G. 



LI SIM 



© — cz: 

Clock time 



sigl — 
PY xor Data 



sig2 > 
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sig3 > 
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Note: Run init.m file first to set parameters 
and generate sigl , sig2, sig3 



Figure 41 LI SIM. M 
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H. L2 SEM.M 
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sigl — 
PY xor Data 
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sig3 — 
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and generate sigl , sig2, sig3 

Figure 42 L2JSIM.M 
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I. 



MPATHSIM.M 



©— 

Clock 

H 

Random 

Number 



t 

time 




x 





"7 




sm 


— w 


jL 

2 



sm 



MATLAB 


| 


Function 


1 



get_delay.m 



-fr [x_delay 
x_delay 



Figure 43 MPATHSIM.M 



93 



J. 



PLL.M 




pll.m 

Figure 44 PLL.M 
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K. 



PLL 2A.M 




Figure 45 PLL2A.M 
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APPENDIX C EXAMPLE SNR CALCULATION 



X 






dB2R(x) = 10 10 convert dB to ratio 


dB(x) =10^1og(x) 


Minimum Received Power 


C:=- 160.0 


~ dB_W 


LI with C/A Code 






Boltzman's Constant 


k =-228.6 


~ dB W/Hz-K 


Noise Temperature 


Tn =28 


- dB K 


Received Noise 


No =k+ Tn 






No =-200.6 


~ dB W/Hz 


Carrier to Noise Density Ratio 


CNR =C- No 






CNR = 40.6 


~ db Hz 



Chip Data Rate 

Information Data Rate 
Band Width - QPSK 

Signal to Noise Ratio 

Processing Gain 
Signal to Noise Ratio 



R = 1.023 10 6 

R_db :=dB(R) 
R_db =60.099 
R_i =50 

B = 0.6R 
B_dB =dB(B) 
B_dB = 57.88 

SNR := CNR- B 
SNR = -17.28 

PG:=dB — ) 

\R_i/ 

SNR:= SNR+ PG 



~ chips/sec 

- dB cps 
~ bits/sec 

~ dB bps 
dB 

~ dB 

PG=43.109 - dB 

SNR = 25.829 



convert ratio to dB 



-dB 
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